% surf_curv_thres.m
% Analyze surface curvature using a threshold.
% The parametric representation of a surface is s = s(u,v)
% where s is a vector function of the two parameters u and v.
% In general each component is a function of u and v so the form is:
% s = sx(u,v)ex + sy(u,v)ey + sz(u,v)ez
% Here we adopt the special case of a Monge patch where
% sx = u, sy = v, and sz = sz(u,v), so the analytical form is: 
% s = uex + vey + f(u,v)ez
% For example, given a field data set of position vectors such as
% (Easting, Northing, Elevation)
% we associate these coordinates with u, v, and f(u,v) respectively.
% This code calculates and plots:
%   surface s in 3D space (x,y,z),
%   principal normal curvature magnitudes, k1 and k2,
%   principal normal curvature directions as vectors on tangent planes,
%   geologic curvature, GC, in 9 different second order classes on s.
% Also calculated are:
%   coefficients of the first fundamental form: E, F, G,
%   coefficients of the second fundamental form: L, M, N,
%   shape operator: I\II.
% Sign conventions:
%   right-hand rule for cross product of tangent vectors (tu x tv) 
%       directs the unit normal vector away from the surface s;
%   curves that are concave upward (pos z) have positive curvature.
% This code was written by David D. Pollard in June 2007 while he was
% at the National Cello Institute (Pomona College).
% Many ideas and snippets of code herein came from earlier codes
%  by Stephan Bergbauer, Ian Mynatt, Steve Martel
% Pete Lovely has contributed corrections to bugs in the June 2007 code. 
% These all are gratefully acknowledged.

clear all; clf % initialize memory and figure

% Read x, y, and z position data for surface from files.
% x and y are horizontal coordinates in a rectangular grid.
% z is vertical coordinate that varies with x and y.
X=dlmread('smax.txt'); Y=dlmread('smay.txt'); Z=dlmread('smaz.txt');
dx=X(1,2)-X(1,1); dy=Y(2,1)-Y(1,1); % Spacing between grid points. 
[m,n]=size(Z); % Number of rows and cols for grid on parameter plane.
surf(X,Y,Z) % Plot wire-frame of s with color fill.
xlabel('x-axis'), ylabel('y-axis'), zlabel('z-axis')
title('raw data surface S=S(u,v)'), axis equal, view([-83,11])

% Begin spectral filtering.
sf=0; % Filter if sf=1; do not filter if sf=0.
if sf==1
    % Spectral filtering removes high frequency content of data
    % using Fast Fourier Transform (fft) methods (reference...).
    % For a digital fft the spatial frequencies are stored as follows:
    % [0, +1df, +2df,..., +(nf-1df), nf, -(nf-1df), ... , -2df, -1df]
    % The spatial frequency axis is discretized into intervals df such that
    % df = 1/(n*ds) where n = # of samples and ds = sample spacing.
    % The minimum non-zero absolute value of the spatial frequency is df.
    % The maximum value is the Nyquist frequency = nf = 1/(2*ds).
    % A spatial frequency threshold, sft, is set such that 
    % all spatial frequencies greater than sft are zeroed out.
    
    sft=0.001;	% Define the spatial frequency threshold.
    % sft should be in the range df < sft < nf
    % If this range is unknown, set sft = 1/eps, run the m-script,
    % identify the range, and set sft within it for subsequent runs.
    % As sft decreases from nf toward df, the surface is smoothed.
    
    Z2=[Z fliplr(Z)]; % Add padding to Z array to reduce edge effects.
    Z3=[Z2;flipud(Z2)];
    [o,p]=size(Z3); % Number of rows and cols for padded Z array.
    dfx=1/(p*dx); % Define interval for the x spatial frequency axis.
    nfx=1/(2*dx); % Define Nyquist frequency for the x axis.
    % Define vector of points on the x spatial frequency axis
    Fx=(1:(p/2))*dfx; Fx=[0, Fx, -fliplr(Fx(1:end-1))];
    dfy=1/(o*dy); % Define interval for the x spatial frequency axis.
    nfy=1/(2*dy); % Define Nyquist frequency for the x frequency axis.
    % Define vector of points on the y spatial frequency axis
    Fy=(1:(o/2))*dfy; Fy=[0, Fy, -fliplr(Fy(1:end-1))];
    TZ=fft2(Z3); % 2D Fourier transform of the padded Z array.
    [FX,FY]=meshgrid(Fx,Fy); % Define 2D grid of spatial frequencies.
    FM=sqrt(FX.^2+FY.^2)<sft; % Define mask, FM, for grid such that:
       % =1 if spatial frequency is < sft
       % =0 if spatial frequency is > sft
    Z4=real(ifft2(FM.*TZ)); % Use mask and inverse fft
    Z=Z4(1:m,1:n); % Select subset of Z array to discard padding.
    figure, surf(X,Y,Z) % Plot wire-frame of s with color fill.
    xlabel('x-axis'), ylabel('y-axis'), zlabel('z-axis')
    title('filtered data surface S=S(u,v)'), axis equal, view([-83,11])
    clear Z2 Z3 Fx Fy TZ FX FY FM Z4
end
% End spectral filtering.

% Begin curvature analysis.
kt = 0.0008; % Define the curvature threshold.
    % In natural data sets principal curvatures are unlikely to be zero.
    % To admit geologic curvature classes in which either the Gaussian
    % or mean curvatures are exactly zero one or both principal
    % curvatures must be zero: principal curvatures are set to zero if
    % their absolute values are less than kt.
    % If the range of absolute values is unknown, set kt = 0, run the
    % m-script, identify the range, and set kt within it for later runs.
    % The value of kt should be based on an analysis of data errors.

SX = X; SY = Y; SZ = Z; % Associate components of s with data values.
% Calculate the first partial derivatives of s w/re to u and v.
% These are components of the tangent vectors to the u- and v-parameter
% curves on s and may be used to plot these vectors.
% Note: for the Monge patch some components are identically 1 or 0.
du=dx; dv=dy; % Distance between grid points on the parameter plane.
SXU=ones(size(SX)); SXV=zeros(size(SX)); % partial derivatives of SX
SYU=zeros(size(SY)); SYV=ones(size(SY)); % partial derivatives of SY
[SZU,SZV]=gradient(SZ,du,dv); % partial derivatives of SZ

% Construct two 3D arrays, each of size (m,n,3).
% Each page in 3rd dimension contains a component of the tangent vectors
% SU and SV are ordered as SXU, SYU, SZU and SXV, SYV, SZV respectively.
SU(:,:,1)=SXU; SU(:,:,2)=SYU; SU(:,:,3)=SZU;
SV(:,:,1)=SXV; SV(:,:,2)=SYV; SV(:,:,3)=SZV;
% Calculate the coefficients of the first fundamental form.
E=dot(SU,SU,3); F=dot(SU,SV,3); G=dot(SV,SV,3);

% Calculate components of the unit normal vector, N, to the surface s.
% Numerator is vector (cross) product of tangent vectors.
% Denominator is absolute value of this cross product.
CUV = cross(SU,SV,3);
AC = sqrt(CUV(:,:,1).^2 + CUV(:,:,2).^2 + CUV(:,:,3).^2);
NX = CUV(:,:,1)./AC; NY = CUV(:,:,2)./AC; NZ = CUV(:,:,3)./AC; 

% Calculate the first partial derivatives of N w/re to u and v.
% These are used to calculate the second fundamental form coefficients.
[NXU,NXV]=gradient(NX,du,dv); [NYU,NYV]=gradient(NY,du,dv); 
[NZU,NZV]=gradient(NZ,du,dv); 

% Construct two 3D arays each of size (m,n,3).
% Each page in 3rd dimension contains a component of the normal vectors.
NU(:,:,1)=NXU; NU(:,:,2)=NYU; NU(:,:,3)=NZU;
NV(:,:,1)=NXV; NV(:,:,2)=NYV; NV(:,:,3)=NZV;
% Calculate the coefficients of the second fundamental form.
L=-dot(NU,SU,3); M=-0.5*(dot(NU,SV,3)+dot(NV,SU,3)); N=-dot(NV,SV,3);

% Preallocate arrays that are built inside the following loops.
K1=zeros(size(SZ)); K1U=zeros(size(SZ)); K1V=zeros(size(SZ));
K2=zeros(size(SZ)); K2U=zeros(size(SZ)); K2V=zeros(size(SZ));

% The first and second fundamental form coefficients are used to
% construct two 2 by 2 matrices, I and II, which then are used to
% calculate the 2 by 2 shape operator matrix, SO, which then is used to
% calculate the principal curvature magnitudes and apparent directions.
for i=1:m % Loop over rows of grid points on parameter plane.
    for j=1:n % Loop over columns of grid points on parameter plane.
        I=[E(i,j) F(i,j);F(i,j) G(i,j)]; % Matrix of first fund. form.
        II=[L(i,j) M(i,j);M(i,j) N(i,j)]; % Matrix of second fund. form.
        SO=inv(I)*II; % 2 by 2 shape operator matrix: same as I\II.
        A=max(max(isnan(SO))); % Test for nan elements in SO.
        if A==1 % if any are nan set all elements of KD and KM to nan
            KD=[NaN NaN;NaN NaN]; KM=[NaN NaN;NaN NaN];
        elseif A==0 % if not nans compute eigenvalues and vectors of SO
            [KD,KM]=eig(SO); % KD and KM are 2 by 2 matrices.
            % Cols of KD are eigenvectors; trace of KM are eigenvalues.
            if KM(1,1)<KM(2,2) % Order principal curvatures so K1 > K2.
                KM=rot90(KM,2); KD=fliplr(KD);
            end
            % KD(:,1) and KM(1,1) are eigenvector and value for K1.
            % KD(:,2) and KM(2,2) are eigenvector and value for K2.
        end
        % Extract principal curvature magnitudes, K1 and K2.
        K1(i,j)=KM(1,1); K2(i,j)=KM(2,2);
        % Extract components of apparent principal curvature directions.
        % Note: these are 2D vectors in the parameter plane which must be
        % projected onto planes tangent to s for actual directions.
        K1U(i,j)=KD(1,1); K1V(i,j)=KD(2,1);
        K2U(i,j)=KD(1,2); K2V(i,j)=KD(2,2);
    end
end

% Project apparent principal curvature vectors onto planes tangent to s
% where they become 3D vectors with components in x, y, and z.
% For the Monge patch some partial derivatives are identically 1 or 0.
K1X=K1U; K1Y=K1V; K1Z=K1U.*SZU+K1V.*SZV;
K2X=K2U; K2Y=K2V; K2Z=K2U.*SZU+K2V.*SZV;
% Normalize direction vector components by vector magnitudes.
K1M = sqrt(K1X.^2+K1Y.^2+K1Z.^2);
K1X = K1X./K1M; K1Y = K1Y./K1M; K1Z = K1Z./K1M;
K2M = sqrt(K2X.^2+K2Y.^2+K2Z.^2);
K2X = K2X./K2M; K2Y = K2Y./K2M; K2Z = K2Z./K2M;

% Plot s with color fill scaled to K1 and superimpose directions.
figure, surf(SX,SY,SZ,K1), colorbar, axis equal, hold on
delz = (max(max(SZ))-min(min(SZ)))/100; % Offset for direction lines.
quiver3(SX,SY,SZ+delz,K1X,K1Y,K1Z,'k.')
title('maximum principal curvature k1 magnitude and direction on s')
xlabel('x-axis'), ylabel('y-axis'), zlabel('z-axis')
% Plot s with color fill scaled to K2 and superimpose directions.
figure, surf(SX,SY,SZ,K2), colorbar, axis equal, hold on
quiver3(SX,SY,SZ+delz,K2X,K2Y,K2Z,'k.');
title('minimum principal curvature k2 magnitude and direction on s')
xlabel('x-axis'), ylabel('y-axis'), zlabel('z-axis')

% Set K1 to 0 if |K1| < kt and set K2 to 0 if |K2| < kt. 
K1(find(abs(K1)<=kt))=0; K2(find(abs(K2)<=kt))=0;
% Calculate Gaussian curvature, KG, and mean curvature, KM.
KG=K1.*K2; KM=0.5*(K1+K2);
% Initialize geological curvature array, GC, to NaN
GC=nan(size(KG));

% Calculate geological curvature, GC, to classify shapes at grid points
% on s into 8 different catagories using the curvature threshold.
% For each catagory set GC equal to the designated color value.
kt2 = kt*kt;
in=find(KG<-kt2 & KM<-kt); % Identify antiformal saddle points.
GC(in)=-1*ones(size(in)); 
in=find(KG==0 & KM<-kt/2); % Identify cylindrical antiformal points.
GC(in)=-2*ones(size(in));
in=find(KG>kt2 & KM<-kt); % Identify domal points.
GC(in)=-3*ones(size(in));
in=find(KG<-kt2 & abs(KM)<=kt); % Identify perfect saddle points.
GC(in)=zeros(size(in)); 
in=find(KG==0 & abs(KM)<=kt/2); % Identify planar points.
GC(in)=4*ones(size(in));
in=find(KG<-kt2 & KM>kt); % Identify synformal saddle points.
GC(in)=1*ones(size(in));
in=find(KG==0 & KM>kt/2); % Identify cylindrical synformal points.
GC(in)=2*ones(size(in));
in=find(KG>kt2 & KM>kt); % Identify basinal points.
GC(in)=3*ones(size(in));

% Plot s with fill color according to geologic curvature table
figure, surf(SX,SY,SZ,GC), axis equal, view([-83,11])
caxis([-4 4]); % Linear color scaling: parameters are cmin and cmax.
colormap(jet); % Choose colormap: default is jet
title('Geological curvature classification of points on s')
xlabel('x-axis'), ylabel('y-axis'), zlabel('z-axis')

% Plot color table for geologic curvature.
% If the parameters for the table are changed above, change here also. 
CX = [1 2 3 4]; CY = [1 2 3 4];
[CXX CYY] = meshgrid(CX,CY); CZZ = zeros(size(CXX));
GCT = [-1 -2 -3; 0 4 nan; 1 2 3];
figure, surf(CXX,CYY,CZZ,GCT), caxis([-4 4]), axis equal
view(0,90), colormap(jet), axis off
text(1.05,1.5,'antiformal saddle            antiform                 dome') 
text(1.15,2.5,'perfect saddle              plane               undefined') 
text(1.05,3.5,'synformal saddle           synform                   basin')
text(1.3,0.85,'KG < -kt^2               KG = 0               KG > kt^2')
text(0.4,1.5,'KM < -kt'), text(0.4,2.5,'|KM| <= kt'), text(0.4,3.5,'KM > kt')
title('Geological curvature with threshold classification color table')
